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Glaucoma is a major cause of blindness in the world. To date, common genetic variants associated with glau- 
coma only explain a small proportion of its heritability. We performed a genome-wide association study of 
intra-ocular pressure (IOP), an underlying endophenotype for glaucoma. The discovery phase of the study 
was carried out in the TwinsUK cohort (N = 2774) analyzing association between IOP and single nucleotide poly- 
morphisms (SNPs) imputed to HapMap2. The results were validated in 12 independent replication cohorts of 
European ancestry (combined N = 22 789) that were a part of the International Glaucoma Genetics 
Consortium. Expression quantitative trait locus (eQTL) analyses of the significantly associated SNPs were per- 
formed using data from the Multiple Tissue Human Expression Resource (MuTHER) Study. In the TwinsUK 
cohort, IOP was significantly associated with a number of SNPs at 9q33.3 (P = 3.48 x 1 0 8 for rs2286885, the 
most significantly associated SNP at this locus), within the genomic sequence of the FAM125B gene. 
Independent replication in a composite panel of 1 2 cohorts revealed consistent direction of effect and significant 
association (P = 0.003, for fixed-effect meta-analysis). Suggestive evidence for an eQTL effect of rs2286885 was 
observed for one of the probes targeting the coding region of the FAM125B gene. This gene codes for a compo- 
nent of a membrane complex involved in vesicular trafficking process, a function similar to that of the Caveolin 
genes (CAV1 and CAV2) which have previously been associated with primary open-angle glaucoma. This study 
suggests a novel association between SNPs in FAM125B and IOP in the TwinsUK cohort, though further studies 
to elucidate the functional role of this gene in glaucoma are necessary. 



INTRODUCTION 

Glaucoma represents a group of optic neuropathies character- 
ized by a typical pattern of optic nerve damage resulting from 
loss of retinal ganglion cells and axons. Glaucoma accounts 
for a significant portion of the global burden of visual impair- 
ment, being the third leading cause of visual impairment and 
the second leading cause of blindness in the world (1). 



Raised intra-ocular pressure (IOP) is a major modifiable risk 
factor for glaucoma (2,3) and treatment directed at lowering 
IOP remains the mainstay of current glaucoma treatment (3,4). 
IOP is a heritable trait, with heritability estimates ranging from 
0 . 3 5 to 0 . 74 in different studies (5 - 7) . Investigating genetic vari- 
ation influencing IOP might help better understand IOP regula- 
tion and its role in glaucoma pathophysiology. So far, in the 



* To whom correspondence should be addressed at: Department of Twin Research and Genetic Epidemiology, St Thomas' Hospital, London SE 1 7EH, UK. 
Tel: +44 2071889055; Email: chris.hammond@kcl.ac.uk 

^See Appendix for the list of International Glaucoma Genetics Consortium members. 
© The Author 2014. Published by Oxford University Press. 

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.Org/licenses/by/3.0/), 
which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. 



3344 Human Molecular Genetics, 2014, Vol. 23, No. 12 



case of glaucoma, case-control studies and quantitative trait- 
based approaches have jointly identified a number of genes 
that influence the risk of developing clinical glaucoma as well 
as cause incremental changes in population-wide risk of the 
underlying quantitative traits. These include the association of 
TMCOl with open-angle glaucoma (OAG) (8) and IOP (9), 
and the association of CDKN2B with OAG (8,10) and vertical 
cup-to-disc ratio (11). This convergence of evidence between 
variants associated with glaucoma and those associated with 
'healthy' variation in quantitative traits underlying the disease, 
justifies the use of quantitative traits (endophenotypes) to 
dissect the genetic basis of glaucoma. 

Here, we describe a genome-wide association study (GWAS) 
investigating the healthy variation of IOP in general populations 
of European ancestry and report findings about genetic variants 
that we found consistently associated with IOP in the participat- 
ing cohorts. 

RESULTS 

Association testing in the TwinsUK cohort was performed for 
~ 1.87 million SNPs that passed our quality control (QC) criteria. 
Both the genomic inflation factor (1.036) and the quantile plot 
(Supplementary material, Fig. SI) of the results indicate no sig- 
nificant inflation of test statistics due to population stratification. 

The region most significantly associated with IOP was 9q3 3 . 3 , 
where three SNPs crossed the conventional threshold for 
genome-wide significance (Figs 1 and 2). In no other region 
SNPs crossed this significance threshold. Of the two previously 
reported IOP loci by van Koolwijk et al. (9), in which TwinsUK 
was one of the participating studies, rsl 1656696 on the GAS7 
locus showed statistically significant replication (P = 1.36 x 
10~ 2 ), although not at a GWAS level; while the association for 
rs7555523 on the TMCOl locus was not statistically significant 
(P= 2.4 x 10 _1 ) (9). 



The most significantly associated SNP at 9q33.3 was 
rs2286885 (P = 3.48 x 10" 8 ). Each copy of the A allele of 
this SNP increased IOP on average by 0.56 mmHg (95% CI: 
0.36-0.75 mmHg). Two other SNPs at this locus that were sig- 
nificantly associated with IOP (rs2286886 and rsl0819169) 
were in perfect linkage disequilibrium (LD) with rs2286885. 
This association remained the same (/3 = 0.56, SE = 0.103, 
P = 6.5 x 10 _8 forrs2286885)afterinclusionofthegenotyping 
batch as a covariate in the analysis. For rs2286885, adjustment 
for central corneal thickness (CCT) improved the evidence for 
its association (/3 on adjustment for CCT = 0.59 mmHg; P = 
1.16 x 10 9 ); adjustment for SBP, on the other hand, decreased 
the significance of association (/3 on adjustment for SBP = 
0.50 mmHg; P = 6.46 x 10~ 7 ). The direction and magnitude 
of effect for rs2286885 remained unchanged on stratifying by 
sex and after selecting one twin from each pair. 

The rs2286885 polymorphism is located within intron 8 of 
the FAM125B gene {Homo sapiens family with sequence 




Chromosome 



Figure 1 . Results of the GWAS in the TwinsUK cohort (The upper line in the plot 
demarcates SNPs that are genome-wide significant (P < 5 x 10 8 ) and the 
lower line demarcates SNPs that show suggestive significance (P < 5 x 1 0 -6 )). 




Figure 2. Regional association plot of the region 9q33.3 in the Twins UK cohort. 
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Table 1. Results for rs2286885 in the discovery and replication studies 



Cohort Number of subjects" Frequency b /3 b SE P-value 



TwinsUK (discovery study) 


2772 


0.55 


0.56 


0.101 


3.48 x 10~ 8 


BATS 


1152 


0.53 


0.147 


0.13 


0.2602 


BMES 


1667 


0.55 


0.0384 


0.099 


0.6994 


ERF 


2589 


0.57 


0.064 


0.089 


0.476 


Framingham 


2455 


0.54 


0.090 


0.091 


0.327 


GHS1 


2727 


0.56 


0.136 


0.078 


0.081 


GHS2 


1130 


0.55 


-0.020 


0.113 


0.858 


ORCADES 


474 


0.54 


0.358 


0.184 


0.051 


RSIII 


2034 


0.57 


0.057 


0.092 


0.537 


RSII 


2116 


0.56 


0.066 


0.094 


0.483 


RSI 


5782 


0.57 


0.101 


0.062 


0.106 


Southampton 


166 


0.55 


0.042 


0.492 


0.932 


TEST 


663 


0.53 


0.010 


0.194 


0.960 



BATS, Brisbane Adolescent Twin Study; BMES, Blue Mountains Eye Study; ERF, Erasmus Rucphen Family Study; GHS, Guttenberg Health Study; ORCADES, 
Orkney Complex Disease Study; RS, Rotterdam Study; TEST, Tasmanian Eye Study of Twins; SE, dtandard error. 
a Number of subjects refers to those with genotype information for rs2286885. 

b Frequency and beta refer to the frequency and the effect size respectively for the A allele of rs2286885. 



less in the replication panel, with each copy of the risk allele 
(A) of rs2286885 increasing IOP by 0.09 mmHg (95% CI: 
0.03-0.14 mmHg). A joint meta-analysis of the discovery and 
replication samples thus reduced the strength of association 
(P = 5.67 x 10" 6 ), with each copy of the risk allele (A) of 
rs2286885 increasing IOP by 0.12 mmHg (95% CI: 0.07- 
0.18 mmHg). 

We checked for associations between rs2286885 and the 
three probes of the HumanHT-12 array covering FAM125B 
coding regions. There was some evidence for an eQTL associ- 
ation, with the best observed SNP-expression association 
being between rs2286885 (which was the top IOP-associated 
SNP at the locus) and the probe (ILMN_1 652525) targeting 
the isoform 1 of FAM1 2 5B in fat tissue (P = 0.02). The observed 
eQTL association however does not survive multiple testing cor- 
rections for number of probes and number of tissues tested. 



DISCUSSION 

Here we report a novel association between variants located on 
the long arm of chromosome 9 and IOP levels. Our conclusions 
were based on GWAS results of one population of European 
descent, and replication of the variants in a compound panel of 
12 independent populations of the same ancestry. The SNPs at 
9q33 .3 associated with IOP are within the LD block that contains 
the gene FAM125B. The protein product of FAM125B forms a 
component of ESCRT-I, a highly conserved complex (12), that 
is involved in vesicular trafficking process (also known as trans- 
cytosis). It is thought that vesicular transport/cell membrane re- 
modelling pathways in the endothelial cells of the eye regulate 
IOP through active modulation of the formation of giant 
vacuoles and endothelial pores (13), which controls the drainage 
rather than the production of aqueous humour in the eye. Further 
support for the importance of this metabolic pathway in the de- 
termination of IOP is the recent discovery of the formation of 
intracellular vesicles in cultured human trabeculocytes in re- 
sponse to treatment with latanoprost, a drug which is a well 
known first line drop treatment to lower IOP (14). Incidentally, 
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Figure 3. Forest plot for rs2286885 for the replication cohorts. ES, effect size; CI, 
confidence intervals; BATS, Brisbane Adolescent Twin Study; BMES, Blue 
Mountains Eye Study; ERF, Erasmus Rucphen Family Study; GHS, Gutenberg 
Health Study; ORCADES, Orkney Complex Disease Study; RS, Rotterdam 
Study; TEST, Tasmanian Eye Study of Twins. 



similarity 125, member B). FAM125B codes for two different 
isoforms (isoform 1 and isoform 2) that contain 10 and 6 
exons, respectively. The protein encoded by this gene (multive- 
sicular body subunit 12B, also known as MVB12B) is a compo- 
nent of the ESCRT-I (endosomal sorting complex required for 
transport I) complex, a heterotetramer, which mediates the 
sorting of ubiquitinated cargo proteins from the plasma mem- 
brane to the endosomal vesicle. 

We attempted to replicate the most significantly associated 
SNP at 9q33.3, rs2286885, in European cohorts participating 
in the International Glaucoma Genetics Consortium (IGGC). 
With the exception of one cohort, the direction of effect for 
rs2286885 was consistent with the one observed in the 
TwinsUK in all the other replication cohorts (Table 1). In a com- 
bined fixed-effect inverse variance meta-analysis of all the rep- 
lication cohorts, rs2286885 was associated with IOP (P = 
0.003) (Fig. 3). The magnitude of the effect was noticeably 
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a recent study has identified association for IOP with variants in 
another gene (ICA1) that is also involved in vesicular transporta- 
tion/cell membrane remodelling (15). In the absence of a case- 
control panel to test the role of the associated SNPs mFAM125B 
in glaucoma, we hypothesize that FAM125B could potentially 
mediate susceptibility to glaucoma through the vesicular trans- 
portation/cell membrane remodelling pathways, which have 
previously been implicated in glaucoma through the reported as- 
sociation of the CA VI - CA V2 genes (16). 

FAM125B mRNA is highly expressed in the human retina 
(http://biogps.org/). The result of our eQTL analysis was con- 
strained by tissue availability. We observed only modest eQTL 
association for rs2286885 and FAM125B in one of the three 
non-ocular tissues we studied; however the possibility of exist- 
ence of such an effect cannot be ruled out as gene expression is 
known to be a tissue-specific phenomenon (17). 

Our results are important for a number of reasons. First they are 
illustrative of the fact that the path towards elucidating the genetic 
architecture of the complex traits such as IOP is still long and that 
there is still place for G WAS to identify genes and genetically con- 
trolled mechanisms that had previously eluded us. The effect size 
of the variants identified in FAM125B is small, as is the portion of 
the phenotypic variance that can be attributed to them. The find- 
ings are consistent with GWAS findings for most complex dis- 
eases to date (18), which has changed the way we look at 
complex disease genetics from an oligogenic to a plurigenic 
model, with hundreds and potentially thousands of individual var- 
iants contributing to phenotypic variation (19). 

Second, our results offer an interesting insight on the phenom- 
enon of 'winner's curse'. Selection of SNP associations, robust 
to conservative Bonferroni correction, to be taken forward to 
subsequent replication stages, often leads to a bias towards 
selecting SNPs whose effect sizes in the discovery cohort are 
over-estimated, albeit truly positive (20). This way of selecting 
replication SNPs encourages a 'winner's curse' phenomenon. 
This is neither a new nor a harmful phenomenon; it has been 
well documented for some of the major confirmed results for a 
variety of traits in the GWAS era (21). Thus, in the discovery 
phase of our study, the estimated effect size of the associated var- 
iants at 9q33.3 appears to be over-estimated; while the effect 
sizes observed in the replication cohorts might reflect the true 
effect of the variants. 

Third, these results show the wide range of possibilities that 
can explain the connection between DNA sequence variations 
with phenotypic expression. The associated variants in 
FAM125B were located within the intron 8 of the gene. The 
exons that flank this intron code for the isoform 1 of the 
protein product; on the other hand, the flanking exons are 
spliced out of the pre-mRNA when coding for the isoform 
2. Moreover, our eQTL results showed association with only 
the isoform 1 of the gene, but not with the isoform 2. Our 
results could suggest a possible link between alternative spli- 
cing/ isoform expression at this locus and IOP regulation, but 
tissue specificity is likely to remain a difficult problem to 
tackle when studying expression of genes involved in ocular phe- 
notypes. 

Fourth, it is known that CCT influences IOP readings (22,23), 
and with CCT being a highly heritable trait (estimates range from 
0.72 to 0.95) (5,6,24), there exists the possibility of CCT con- 
founding genetic associations for IOP. Adjustment for CCT in 



our analysis improved the strength of association for SNPs at 
9q33.3 locus, thus highlighting the fact that the locus regulates 
IOP rather than CCT. This also supports the notion that pheno- 
type refinement for complex traits offers an effective strategy 
to improve the power to detect genetic variants underlying them. 

A potential drawback of our study is that the discovery phase 
was conducted in a cohort with specific demographic features — 
a UK-based, predominantly female population. This might limit 
generalization, so replication of our findings in other cohorts 
with differing demographic structures is important. Three out 
of the four Asian cohorts of IGGC did not replicate the direction 
of effect for rs2286885 (Supplementary Material, Table S2), a 
finding that could reflect differences in LD pattern at this locus 
between European and Asian populations, justifying our deci- 
sion to include only the European cohorts in the replication. 
We further investigated the population diversity of the asso- 
ciated variants at 9q33.3 in European (HapMap-CEU) and 
Asian (HapMap-HCB and HapMap-JPT) population panels of 
HapMap. A significant difference in the allele frequency 
pattern of rs2286885 was noted between European and Asian 
populations — the frequency of the risk allele (A) of this SNP 
was about 55% in the European population, while it was 25 
and 28% in the two Asian populations. This is indicative of dif- 
fering haplotypic frequencies, and thus differences in the pat- 
terns of LD at this locus between European and Asians 
populations. Such variations in LD patterns between popula- 
tions, if not accounted for, can confound the ability to detect 
associations when analyzing ethnically heterogenous population 
panels (25). 

Although GWAS evidence presented here suggests associ- 
ation between FAM125B and IOP, short of any direct biological 
evidence, association results remain probabilistic that still have 
the potential of a type I error. Results from individual GWAS 
may be opening small windows into the genetic architecture of 
IOP and, by extension, glaucoma; however, functional charac- 
terization of the gene function will be necessary to fulfil the po- 
tential translational benefits of such studies. 



MATERIALS AND METHODS 

We analysed 2774 participants (95% female and all of Caucasian 
ancestry) within the TwinsUK adult twin registry based at 
St. Thomas' Hospital in London (26) for whom both genotype 
and IOP information was available (Supplementary Material, 
Table SI). Of the 2774 subjects, 681 were monozygotic twin 
pairs, 62 7 were dizygotic twin pairs and 158 were unrelated indi- 
viduals. Twins largely volunteered unaware of the eye studies 
interests at the time of enrolment and gave fully informed 
consent under a protocol reviewed by the St. Thomas' Hospital 
Local Research Ethics Committee. Exclusion criteria included 
any form of glaucoma surgery such as trabeculectomy or laser 
surgery that could alter IOP. 

We measured IOP with a non-contact air-puff tonometer. The 
Ocular Response Analyser (ORA; Reichert®, Buffalo, NY, 
USA) ejects an air impulse in order to flatten the cornea, which 
is detected by an electro-optical collimation system. The mean 
IOP was calculated from four readings (two from each eye) for 
each participant. IOP for subjects receiving IOP-lowering med- 
ications (26 out of 2774) was imputed by increasing the 
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measured value by 30%, based on efficacy data from commonly 
prescribed therapies (27). As CCT and systolic blood pressure 
(SBP) are known to influence IOP measurements (23), they 
were included as covariates in the analysis. CCT was measured 
using an ultrasound pachymetry device provided with the ORA 
instrument. Three SBP measurements were made for each 
subject using automated calibrated instrument, of which the 
average of the second and the third readings were considered 
for the analysis. 

Subjects were genotyped in two different batches of approxi- 
mately the same size, using two genotyping platforms from Illu- 
mina: 300 K Duo and HumanHap6 1 0-Quad arrays. Whole 
genome imputation of the genotypes was performed using 
HapMap2 (www.hapmap.org) haplotypes. 

Stringent QC measures were implemented, including 
minimum genotyping success rate (>95%), Hardy-Weinberg 
equilibrium (P > 10 _6 ),minimumMAF(> 1%) and imputation 
quality score (>0.7). Subjects of non-Caucasian ancestry were 
excluded from the analysis. 

A linear regression model, adjusted for age and sex, was fitted 
to test for association between genome-wide SNPs as independ- 
ent and IOP as the dependent variable. An additive model of 
effect for the risk allele of an SNP was implemented. Additional 
analyses were performed with CCT and SBP as covariates in the 
analysis. A score test statistic as implemented in MERLIN (28), 
that takes into account the pedigree structure and the zygosity of 
the twins, was used to adjust for the non-independence of the 
observations. 

Loci conventionally considered genome-wide significant 
(P < 5 x 10~ 8 ) in the discovery cohort were taken forward 
and meta-analysed using the summary statistics data obtained 
from 12 independent cohorts of European ancestry that were 
a part of the IGGC. Individual cohorts performed the replication 
association analyses for rs2286885 based on a protocol decided 
by the consortium in order to ensure consistency in the analyses. 
The replication cohorts had a combined sample size of 22 789; 
complete descriptions of the study populations, phenotyping 
and genotyping methods for the replication studies are provided 
in the Supplementary Material. Since differences in LD patterns 
between populations are known to affect the portability of 
phenotypic associations when the replication effort is attempted 
in populations that are distinct from the original population in 
which the genome-wide study is performed (25), four cohorts 
of Asian ancestry that were also a part of the IGGC were not 
included in our replication analysis. Where more than one SNP 
at a locus was genome-wide significant, i.e. in presence of LD, 
only the single most associated SNP was chosen for replication. 
A fixed-effect inverse variance meta-analysis of all the cohorts 
was performed using the module 'metan' on Stata Statistical 
Software, 1 1 (College Station, TX, USA). 

Gene expression data for a subset of the TwinsUK cohort was 
obtained from the MuTHER study (29). As a part of the 
MuTHER study, 855 subjects from the TwinsUK cohort had 
their transcript expression quantified in three different tissue 
types (skin, fat and LCLs). This was done for 18 170 genes 
across the genome using 27 499 probes on Illumina's whole 
genome expression array (HumanHT-12 version 3) further 
details of the study methods are available in Nica et al. (29). 
For the genome-wide significant SNPs in the TwinsUK cohort, 
we tested for a possible eQTL effect on the overlapping gene. 



SUPPLEMENTARY MATERIAL 

Supplementary Material is available at HMG online. 
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